#ifdef  SEQ_ENDBG
#define DBG(...)		printf( __VA_ARGS__ )
#else
#define DBG(...)
#endif

#include <stdio.h>
#include <assert.h>
#include <omp.h>
#include <alloca.h>

#include "bitonic_interface.h"

#define SWAP(a, b)  \
			a ^= b; \
			b ^= a; \
			a ^= b

#define CMP_SWAP(a, b)     {\
			if(a > b) {     \
				SWAP(a, b); \
			}               \
		}


static void bitonic_merge_rec(data_t*, size_t, size_t);
void bitonic_merge_itr(data_t*, size_t);

/* Debugging */
extern void print_array(data_t*, size_t);

//static
void bitonic_merge_itr(data_t *data, size_t l) {
	DBG("Starting iteration\n");
	if(l <= 1) {
		return;
	}

	uint8_t *odd = alloca(l);
	for(size_t t = 0; t < l; t++) {
		odd[t] = 0;
	}

	size_t tmp = l;
	size_t log_l = 0;
	while(tmp >>= 1) {
		log_l++;
	}
	if(l%2) {
		log_l++;
	}

	DBG("log[%ld] l[%ld]\n", log_l, l);

	for(size_t k = 1; k <= log_l; k++) {
		/* k is the depth of the recursion tree */
		size_t step;
		size_t l_parent = (l/(1 << (k-1)));

		print_array(data, l);
		DBG(" k[%ld] parent[%ld]\n", k, l_parent);

		for(size_t m = 0; m < l; m += l_parent + odd[m]) {
			step = (l/(1 << k));

			DBG("  m[%ld] l_p[%ld]", m, l_parent+odd[m]);
			if((l_parent + odd[m])%2) {
				if(data[ m+step ] < data[ m+step*2 ]) {
					DBG(" odd[%ld]: 1", m);
					odd[m] = 1;
				} else {
					DBG(" odd[%ld]: 1", m+step);
					odd[m+step] = 1;
				}
				step++;
			}
			DBG(" step[%ld]\n", step);

			//DBG("  i = m[%ld] -> [%ld] + [%d]\n", m, m+(l/(1 << k)), odd[m]);
			for(size_t i = m; i < m+(l/(1 << k))+odd[m]; i++) {
				//DBG("   i[%ld] j[%ld]\n", i, i+step);
				CMP_SWAP(data[i], data[i+step]);
			}
		}

		//print_array(data, l);
		//DBG("[0] +(%ld), [%ld] +(%ld)\n", q, q, l-q);
	}

	//bitonic_merge_rec(data, q, 0);
	//bitonic_merge_rec(data+q, l-q, 0);
	DBG("Finishing iteration\n");
}

static
void bitonic_merge_rec(data_t *data, size_t l, size_t n) {
	if(l <= 1) {
		return;
	}

	size_t step = (l/2);
	size_t q = step;

	if(l%2 == 1) {
		if(data[ step ] < data[ (l-1) ]) {
			q++;
		}
		step++;
	}


	//print_array(data, l);
	DBG("[0] +(%ld), [%ld] +(%ld)\n", q, q, l-q);

	if(n > 1) {
		#pragma omp parallel for
		for(size_t i = 0; i < (l/2); i++) {
			CMP_SWAP(data[i], data[i+step]);
		}
		//printf ("parent: id = %d, n = %ld\n", omp_get_thread_num(), n);
		#pragma omp parallel sections num_threads(2)
		{
			#pragma omp section
			{
				//printf ("id = %d, n = %ld\n", omp_get_thread_num(), n/2);
				bitonic_merge_rec(data, q, n/2);
			}
			#pragma omp section
			{
				//printf ("id = %d, n = %ld\n", omp_get_thread_num(), n/2);
				bitonic_merge_rec(data+q, l-q, n/2);
			}
		}
	} else {
		for(size_t i = 0; i < (l/2); i++) {
			CMP_SWAP(data[i], data[i+step]);
		}
		bitonic_merge_rec(data, q, 0);
		bitonic_merge_rec(data+q, l-q, 0);
	}
}

void bitonic_merge(data_t *data, size_t q, size_t l, int p) {
	if(l == 1) {
		return;
	}
	size_t m = (l-1)/2;
	DBG(" q[%ld] l[%ld]\n", q, l);

	/* [0] -> [q], [q+1] -> [l-1] */
	if(q <= m) {
		/*size_t i = 0;
		size_t j = (2*q)+1;

		for(; i <= q; i++, j--) {
			DBG("  0:0  i[%ld] j[%ld]\n", i, j);
			CMP_SWAP(data[i], data[j]);
		}*/

		#pragma omp parallel for
		for(size_t k = 0; k <= q; k++) {
			//DBG("  -:-  i[%ld] j[%ld]\n", k, (((2*q)+1)-k));
			CMP_SWAP(data[k], data[((2*q)+1)-k]);
		}
	} else /* q > m */ {
		/*size_t i = q-((l-1)-q)+1;
		size_t j = (l-1);

		for(; j > q; i++, j--) {
			DBG("  1:1  i[%ld] j[%ld]\n", i, j);
			CMP_SWAP(data[i], data[j]);
		}*/

		/* i = ((l-1)-k)+q-((l-q)-q)+1)
		 *   = 2q-k+1 */
		#pragma omp parallel for
		for(size_t k = (l-1); k > q; k--) {
			//DBG("  -:-  i[%ld] j[%ld]\n", (2*q)-k+1, k);
			CMP_SWAP(data[(2*q)-k+1], data[k]);
		}
	}

	// number of threads
	size_t n;
	if(p == 0) {
		n = omp_get_max_threads();
	} else {
		n = p;
	}

	//bitonic_merge_itr(data, 5);//q+1);
	//bitonic_merge_itr(data+q+1, l-q-1);
	//return;

	omp_set_nested(1);
	//printf(" parent: id = %d, n = %ld\n", omp_get_thread_num(), n);
	#pragma omp parallel sections num_threads(2)
	{
		#pragma omp section
		{
			//printf("  0 id = %d, n = %ld\n", omp_get_thread_num(), n/2);
			bitonic_merge_rec(data, q+1, n/2);
		}
		#pragma omp section
		{
			//printf("  1 id = %d, n = %ld\n", omp_get_thread_num(), n/2);
			bitonic_merge_rec(data+(q+1), l-q-1, n/2);
		}
	}
}
